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ANOVA FOR LONGITUDINAL DATA WITH MISSING VALUES 1 

By Song Xi Chen and Ping-Shou Zhong 

Iowa State University and Peking University, and Iowa State University 

We carry out ANOVA comparisons of multiple treatments for 
longitudinal studies with missing values. The treatment effects are 
modeled semiparametrically via a partially linear regression which 
is flexible in quantifying the time effects of treatments. The empiri- 
cal likelihood is employed to formulate model-robust nonparametric 
ANOVA tests for treatment effects with respect to covariates, the 
nonparametric time-effect functions and interactions between covari- 
ates and time. The proposed tests can be readily modified for a va- 
riety of data and model combinations, that encompasses parametric, 
semiparametric and nonparametric regression models; cross-sectional 
and longitudinal data, and with or without missing values. 

1. Introduction. Randomized clinical trials and observational studies are 
often used to evaluate treatment effects. While the treatment versus control 
studies are popular, multi-treatment comparisons beyond two samples are 
commonly practised in clinical trails and observational studies. In addition to 
evaluate overall treatment effects, investigators are also interested in intra- 
individual changes over time by collecting repeated measurements on each 
individual over time. Although most longitudinal studies are desired to have 
all subjects measured at the same set of time points, such "balanced" data 
may not be available in practice due to missing values. Missing values arise 
when scheduled measurements are not made, which make the data "unbal- 
anced." There is a good body of literature on parametric, nonparametric 
and semiparametric estimation for longitudinal data with or without miss- 
ing values. This includes Liang and Zeger (1986), Laird and Ware (1982), 
Wu, Chiang and Hoover (1998), Wu and Chiang (2000), Fitzmaurice, Laird 
and Ware (2004) for methods developed for longitudinal data without miss- 
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ing values; and Little and Rubin (2002), Little (1995), Laird (2004), Robins, 
Rotnitzky and Zhao (1995) for missing values. 

The aim of this paper is to develop ANOVA tests for multi-treatment 
comparisons in longitudinal studies with or without missing values. Sup- 
pose that at time t, corresponding to k treatments there are k mutually 
independent samples, 

where the response variable Yji(t) and the covariate Xji(t) are supposed to 
be measured at time points t = tjn, . . . ,tjiTj- Here Tj is the fixed number 
of scheduled observations for the jth treatment. However, {Yji(t), Xj^t)} 
may not be observed at some times, resulting in missing values in either the 
response Yji(t) or the covariates Xji(t). 

We consider a semiparametric regression model for the longitudinal data 

(1 x) Yji(t) = XjMPjo + M T (X ji (t),t)j j0 + g j0 (t) + Sji(t), 

j = 1, 2, . . . , k, 

where M(Xji(t),t) are known functions of Xji(t) and time t representing 
interactions between the covariates and the time, /3jo and jjo are p- and q- 
dimensional parameters, respectively, gjo{t) are unknown smooth functions 
representing the time effect, and {sji(t)} are residual time series. Such a 
semiparametric model may be viewed as an extended partially linear model. 
The partially linear model has been used for longitudinal data analysis; see 
Zeger and Diggle (1994), Zhang et al. (1998), Lin and Ying (2001), Wang, 
Carroll and Lin (2005). Wu, Chiang and Hoover (1998) and Wu and Chi- 
ang (2000) proposed estimation and confidence regions for a semiparametric 
varying coefficient regression model. Despite a body of works on estima- 
tion for longitudinal data, analysis of variance for longitudinal data have 
attracted much less attention. A few exceptions include Forcina (1992) who 
proposed an ANOVA test in a fully parametric setting; and Scheike and 
Zhang (1998) who considered a two sample test in a fully nonparametric 
setting. 

In this paper, we propose ANOVA tests for differences among the /3jo' s 
and the baseline time functions <?jo' s ) respectively, in the presence of the 
interactions. The ANOVA statistics are formulated based on the empirical 
likelihood [Owen (1988, 2001)], which can be viewed as a nonparametric 
counterpart of the conventional parametric likelihood. Despite its not re- 
quiring a fully parametric model, the empirical likelihood enjoys two key 
properties of a conventional likelihood, the Wilks' theorem [Owen (1990), 
Qin and Lawless (1994), Fan and Zhang (2004)] and Bartlett correction 
[DiCicco, Hall and Romano (1991), Chen and Cui (2006)]; see Chen and 
Van Keilegom (2009) for an overview on the empirical likelihood for regres- 
sion. This resemblance to the parametric likelihood ratio motivates us to 
consider using empirical likelihood to formulate ANOVA test for longitu- 
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dinal data in nonpar ametric situations. This will introduce a much needed 
model-robustness in the ANOVA testing. 

Empirical likelihood has been used in studies for either missing or lon- 
gitudinal data. Wang and Rao (2002), Wang, Linton and Hardle (2004) 
considered an empirical likelihood inference with a kernel regression impu- 
tation for missing responses. Liang and Qin (2008) treated estimation for 
the partially linear model with missing covariates. For longitudinal data, 
Xue and Zhu (2007a, 2007b) proposed a bias correction method to make 
the empirical likelihood statistic asymptotically pivotal in a one sample par- 
tially linear model; see also You, Chen and Zhou (2006) and Huang, Qin 
and Follmann (2008). 

In this paper, we propose three empirical likelihood based ANOVA tests 
for the equivalence of the treatment effects with respect to (i) the covari- 
ate Xji] (ii) the interactions M(Xji(t),t) and (hi) the time effect functions 
<7jo(0' s ) by formulating empirical likelihood ratio test statistics. It is shown 
that for the proposed ANOVA tests for the covariates effects and the in- 
teractions, the empirical likelihood ratio statistics are asymptotically chi- 
squared distributed, which resembles the conventional ANOVA statistics 
based on parametric likelihood ratios. This is achieved without parametric 
model assumptions for the residuals in the presence of the nonparametric 
time effect functions and missing values. Hence, the empirical likelihood 
ANOVA tests have the needed model-robustness. Another attraction of the 
proposed ANOVA tests is that they encompass a set of ANOVA tests for 
a variety of data and model combinations. Specifically, they imply specific 
ANOVA tests for both cross-sectional and longitudinal data; for parametric, 
semiparametric and nonparametric regression models; and with or without 
missing values. 

The paper is organized as below. In Section 2, we describe the model and 
the missing value mechanism. Section 3 outlines the ANOVA test for com- 
paring treatment effects due to the covariates: whereas the tests regarding 
interaction are proposed in Section 5. Section 4 considers ANOVA test for 
the nonparametric time effects. The bootstrap calibration to the ANOVA 
test on the nonparametric part is outlined in Section 6. Section 7 reports 
simulation results. We applied the proposed ANOVA tests in Section 8 to 
analyze an HIV-CD4 data set. Technical assumptions are presented in the 
Appendix. All the technical proofs to the theorems are reported in a sup- 
plement article [Chen and Zhong (2010)]. 

2. Models, hypotheses and missing values. For the ith individual of the 
jth treatment, the measurements taken at time tji m follow a semiparametric 
model 



(2.1) 
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for j = 1, . . . , k, i = 1, . . . , rij, m = 1, . . . , Tj. Here /3jo and -jjq are unknown 
p- and (/-dimensional parameters and gjo(t) are unknown functions repre- 
senting the time effects of the treatments. The time points {ijim}m=i are 
known design points. For ease of notation, we write (Yji m , Xj im , Mj im ) to 
denote (Y j i(t j i m ),Xj i (t jim ),M T (X ji (t j i m ),t jim )). Also, we will use X T jim = 
(Xj im , MJ im ) and £J = (/3J,tJ). For each individual, the residuals {eji(t)} 
satisfy EfaitylXjift)} = 0, Var{ £ii (t)|X,i(t)} = a?(t) and 

Cov{£j i (t),£ji(s)|X 7 i(t),Xji(s)} = p j (s,t)a j (t)a j (s), 

where pj(s,t) is the conditional correlation coefficient between two residuals 
at two different times. And the residual time series {eji(t)} from different 
subjects and different treatments are independent. Without loss of general- 
ity, we assume t,s £ [0, 1]. For the purpose of identifying f3jQ, 7^0 and gjo(t), 
we assume 

(Pjo,ljO,9jo) = argmin— r ^ ^ E{Y jim - Xj im f3j - Mj im ^ - gj (t jim )} 2 . 

(Pj,H,9i) Tl 3 1 l i=l m=l 

We also require that ^ YaLi EHi s (^imXJ im ) > °> where X jim = X jim - 
E(Kji m \tji m ). This condition also rules out M(Xji(t),t) being a pure func- 
tion of t, and hence it has to be genuine interaction. For the same reason, 
the intercept in model (2.1) is absorbed into the nonparametric part gjo(t). 

As commonly exercised in the partially linear model [Speckman (1988); 
Linton and Nielsen (1995)], there is a secondary model for the covariate Xji m : 

Xji m — hj (tjim) ^jirnj 

(2-2) 

j — 1,2, ... ,k,i — 1, ... , rij, m — 1, . . . ,Tj , 

where /ij(-)' s are p-dimensional smooth functions with continuous second 
derivatives, the residual Uj im = («j im , ■ ■ ■ , u p - im ) T satisfy E(uj im ) = and Uji 
and Ujk are independent for I ^ k, where Uji = (v,ju, . . . , UjiTj)- By the iden- 
tification condition given above, the covariance matrix of Uji m is assumed 
to be finite and positive definite. 

We are interested in testing three ANOVA hypotheses. The first one is on 
the treatment effects with respect to the covariates: 

H 0a : /3io = #20 = • ■ • = Pko vs. H la : f3 i0 / /3 j0 for some i / j. 

The second one is regarding the time effect functions: 

Hob-gw(-) = ■■■ = 9ko{-) vs. H lb :g i0 (-) ^S'jo(-) for some i/j. 

The third one is on the existence of the interaction Hq c : jjq = and H\ c : 
7j0 / 0. And the last one is the ANOVA test for 

H od : 710 = 720 = • • • = 7fc0 vs. H ld : -y i0 / 7 j0 for some i / j. 
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Let Xji = {Xjio,...,Xj iTj } and Y j{ = {Y ji0 , . . . ,Y jiTj } be the complete 
time series of the covariates and responses of the (j,i)th subject (the ith 
subject in the jth treatment), and Y jitid = {Y ji{ t_ d ) , . . . , Y^(t_i) } and X jitA = 
{Xj^u-d)-, ■ ■ ■ be the past d observations at time t for a positive 

integer d < min, {Tj}. For t < d, we set d = t — 1. 

Define the missing value indicator Sju = 1 if (XJ it ,Yjn) is observed and 
Sjit = if (Xj it ,Yjit) is missing. Here, we assume Xjn and Yja are either 
both observed or both missing. This simultaneous missingness of Xja and 
Yjit is for the ease of mathematical exposition. We also assume that 5jio = 1, 
namely the first visit of each subject is always made. 

Monotone missingness is a common assumption in the analysis of lon- 
gitudinal data [Robins, Rotnitzky and Zhao (1995)]. It assumes that if 
Sji(t-i) = then 5ju = 0. However, in practice after missing some scheduled 
appointments people may rejoin the study. This kind of casual drop-out ap- 
pears quite often in empirical studies. To allow more data being included 
in the analysis, we relax the monotone missingness to allow segments of 
consecutive d visits being used. Let <5jit,d = 11?= i $ji(t-i) ■ We assume the 
missingness of (Xj it ,Yjit) is missing at random (MAR) Rubin (1976) given 
its immediate past d complete observations, namely 

P{fijit = M^jit.d = l,Xji,Yji) = P{8jit = M5ji t A = l,Xji*fi,Yji tl i) 

(2-3) _ _ 

= pj {Xj it ^, Yj it ^] Ojo) . 

Here the missing propensity pj is known up to a parameter Ojo. To allow 
derivation of a binary likelihood function, we need to set Sjn = if 6jn t d = 
when there is some drop-outs among the past d visits, which is only tem- 
porarily if 6 jn = 1. This set-up ensures 

(2.4) P(8 ju = 0\S jit>d = 0, X ja>d , Y jit>d ) = 1. 

Now the conditional binary likelihood for {$jit\t=i gi ven Xji an d ^ji is 
i0> • • • j fijiTj Yji) 
Tj 

= J^J P{fijim\fiji(m—1) > ■ ■ ■ > ^jiOi Xji, Yji) 
m=l 

T ' 

— J~J P(&jim \&jim,d — 1? -^jim,dt Yjim,d) 
m=l 

= 1 1 ^Pj (Xjirn,di Yjim,d'i @j ) JIm { 1 — Pj {Xjim,di Yjim,d'i 



m=l 
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In the second equation above, we use both the MAR in (2.3) and (2.4). 
Hence, the parameters Ojo can be estimated by maximizing the binary like- 
lihood 



CniO,) llll/M-V;-/,/.^..,/:^ 



i=lt=l 

(2-5) ^ _ 

x {1 -pjiXj^Yj^O^ 1 - 8 ^]***'*. 

Under some regular conditions, the binary maximum likelihood estimator 
9j is -y/n-consistent estimator of 9jo; see Chen, Leung and Qin (2008) for 
results on a related situation. Some guidelines on how to choose models for 
the missing propensity are given in Section 8 in the context of the empiri- 
cal study. The robustness of the ANOVA tests with respect to the missing 
propensity model are discussed in Sections 3 and 4. 

3. ANOVA test for covariate effects. We consider testing for Ho a : /?io = 
P20 = ■ ■ ■ = Pko with respect to the covariates. Let ir jim (9j) = ]XlL m -dPA X jtt,di 
Yjil,d',6j) be the overall missing propensity for the (j,i)th subject up to 
time tji m . To remove the nonparametric part in (2.1), we first estimate the 
nonparametric function gjo(t). If (3jo and jjo were known, gjo(t) would be 
estimated by 



(3.1) gj(t; p jQ ) = 2^ 2^ w jim , h (t){Y jim - Xj im /3 j0 - MJ im -y j0 ), 

i=l m=l 

where 

(3.2) Wjim,hj{t) = —^—^r. — 



YZL uUSjsi/^siie^Kh. (t jsl - 1) 

is a kernel weight that has been inversely weighted by the propensity KjimiPj) 
to correct for selection bias due to the missing values. In (3.2), K is a uni- 
variate kernel function which is a symmetric probability density, (t) = 
K(t/hj)/hj and hj is a smoothing bandwidth. The conventional kernel es- 
timation of gjo(t) without weighting by 7Tj s i(&j) may be inconsistent if the 
missingness depends on the responses Yju , which can be the case for missing 
covariates. 

Let Aji m denote any of Xji m ,Yji m and Mji m and define 

nj Tj 

(3.3) Aji m — Aji m ^ ^ ^ ^ iVji irnit i l j(tji m s )Aji imi 

ii=l mi=l 
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to be the centering of A^ m by the kernel conditional mean estimate, as is 
commonly exercised in the partially linear regression [Hardle, Liang and Gao 
(2000)]. An estimating function for the (j,i)th subject is 



Tj 6 



m=l ' K 3 im 



where jj is the solution of 



3 ^3 c 

E E —^M jim (Y jim - X] m (3 j0 - MJimJj) = 



i=\ m =l Kjirnyvj 



at the true /3jQ. Note that E{Zji((3jo)} =o(l). Although it is not exactly 
zero, Zji((3jo) can still be used as an approximate zero mean estimating 
function to formulate an empirical likelihood for /?,• as follows. 

Let {pji}™=i be nonnegative weights allocated to {(Xj^ Y^)}^. The em- 
pirical likelihood for f3j is 

(3.4) L nj (/3 i ) = max|f| Pii J, 

subject to Y^LiPji = 1 and YZiLxVjiZjiiPj) =0. 

By introducing a Lagrange multiplier Xj to solve the above optimization 
problem and following the standard derivation in empirical likelihood [Owen 
(1990)], it can be shown that 

(3 - 5) ^fcir^W* 



where Aj satisfies 



(3.6) 



\ ^ Z ji(Pj) 

+ 



The maximum of L nj (/3j) is n^ii n 7 ' achieved at /3j = f3j and Xj = 0, where 
solves YaLi z jiWj) = °- 

Let n = ^2i = \ n j, n j / n Pj for some nonzero pj as n — > oo such that 
X^i=i Pj = 1- As the k samples are independent, the joint empirical likelihood 
for (Pi,P 2 , ■■■,Pk) is 

k 

L n (p 1 ,{3 2 ,...,{3 k ) = l[L n] ({3 j ). 

j'=i 
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The log likelihood ratio statistic for Ho a is 

k 

£ n := -2maxlogL n ,(/3,/3,. . . ,/3) + ^njlogrij 

13 3=1 

(3.7) 

k rij 

= 2min^^log{l + AJZ Jl (/3)}. 

j=l i=l 

Using a Taylor expansion and the Lagrange multiplier to carry out the 
minimization in (3.7), the optimal solution to f3 is 

(3.8) (E^Vn"*) (XX 5 7X^ + °^)' 

where Bj = lim^-^ (n/Z))~ Yn=i E { z ji(Pjo) z jiWjo) T }, 

1 nj Tj f 5 



and 



1 <5- 

= /TTr" E E ^ /a ,Xji m (Yjim - MJimfJj). 
\J n 3 L d i= l m =l KjimiVj) 

The ANOVA test statistic (3.7) can be viewed as a nonparametric coun- 
terpart of the conventional parametric likelihood ratio ANOVA test statistic, 
for instance that considered in Forcina (1992). Like its parametric counter- 
part, the Wilks theorem is maintained for £ n . 

Theorem 1. If conditions A1-A4 given in the Appendix hold, then un- 
der H 0a , 



The theorem suggests an empirical likelihood ANOVA test that rejects 
H 0a if t n > X(jfc_i) PjQ! where a is the significant level and X( fe _ 1)p Q , is the 
upper a quantile of the X?fe_i) p distribution. 

We next evaluate the power of the empirical likelihood ANOVA test under 
a series of local alternative hypotheses: 

— 1/2 

Hi a : Pjo = Pio + c n n j for 2 < j < k, 

where {c n } is a sequence of bounded constants. Define = (/3J" — /3J , j3\ — 
• • • 5 /3i r o " /5fco) r ; D H = tt^iw " '-K^-h,,,, for 2 < j < k and D = 
{D T 12 ,Dl z ,... 1 D T lk ) T . Let E D =Var(D) and 7 2 = A^ET^A^. Theorem 2 
gives the asymptotic distribution of £ n under the local alternatives. 
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Theorem 2. Suppose conditions A1-A4 in the Appendix hold, then un- 
der Hi a , i n A X( fc _ 1)p (7 2 ) as n^oo. 

It can be shown that 

(3.9) s D = n-^B^i^ ® i (fc _ 1} + dm g {n^B 2 n-^ . . . , o-^feO- 1 }. 

As each fi" 1 is (^(n 1 / 2 ), the noncentral component 7 2 is nonzero and bounded 
The power of the a level empirical likelihood ANOVA test is 

0(7) = - p {X( fc -i)p(7 2 ) > X 2 ( k -i)p,J- 

This indicates that the test is able to detect local departures of size 0(re~ 1 / 2 ) 
from Ho a , which is the best rate we can achieve under the local alternative 
set-up. This is attained despite the fact that nonpar ametric kernel estima- 
tion is involved in the formulation, which has a slower rate of convergence 
than y/n, as the centering in (3.3) essentially eliminates the effects of the 
nonparametric estimation. 

Remark 1. When there is no missing values, namely all 5ji m = 1, we 
will assign all 7Tji m (8j) = 1 and there is no need to estimate each OjQ. In this 
case, Theorems 1 and 2 remain valid. It is a different matter for estimation 
as estimation efficiency with missing values will be less than that without 
missing values. 

Remark 2. The above ANOVA test is robust against misspecifying the 
missing propensity Pj(-;9jo) provided the missingness does not depend on 
the responses Yjit,d- This is because despite the mispecification, the mean of 
Zji(f3) is still approximately zero and the empirical likelihood formulation 
remains valid, as well as Theorems 1 and 2. However, if the missingness 
depends on the responses and if the model is misspecified, Theorems 1 and 2 
will be affected. 

Remark 3. The empirical likelihood test can be readily modified for 
ANOVA testing on pure parametric regressions with some parametric time 
effects gjo(t;r]j) with parameters rjj. When there is absence of interaction, 
we may formulate the empirical likelihood for (f3j,rjj) £ R p+S using 

m=1 TTjimWj) \ OT lj J 

x {Yjim ~ Xji m (3j — gjoitjirniVj)} 

as the estimating function for the (j, i)th subject. The ANOVA test can 
be formulated following the same procedures from (3.5) to (3.7), and both 
Theorems 1 and 2 remaining valid after updating p with p + s where s is 
the dimension of rjj . 
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In our formulation for the ANOVA test here and in the next section, 
we rely on the Nadaraya- Watson type kernel estimator. The local linear 
kernel estimator may be employed when the boundary bias may be an issue. 
However, as we are interested in ANOVA tests instead of estimation, the 
boundary bias does not have a leading order effect. 

4. ANOVA test for time effects. In this section, we consider the ANOVA 
test for the nonparametric part 

Hob-gio(-) = ■■■ =9ko(-)- 

We will first formulate an empirical likelihood for gjo(t) at each t, which 
then lead to an overall likelihood ratio for i?ofc- We need an estimator of 
gjo(t) that is less biased than the one in (3.1). Recall the notation defined in 
Section 2: XJ im = (XJ im ,MJ im ) and £J = (/3J,7j). Plugging-in the estimator 

£j to (3.1), we have 

( 4 - 1 ) Jj(*)=J]J] w jimAj (t) (Y jim - X T jim ij ) . 

i=l m=l 

It follows that, for any t £ [0, 1], 

9j{t) ~ gjo(t) = 22 W jiin,h 3 {-t){£ji{tjirn) + - £j) 

i=l m=l 

(4.2) 

+ gjo{tjim) — 9jo(t)}- 

However, there is a bias of order in the kernel estimation since 

Yl ^2 w jim,hj(t){gjo(tjim) - gjo(t)} = \^J z 2 K (z) dz^g'-^t)^ + o p {h^) . 

If we formulated the empirical likelihood based on (jj(t), the bias will con- 
tribute to the asymptotic distribution of the ANOVA test statistic. To avoid 
that, we use the bias-correction method proposed in Xue and Zhu (2007a) 
so that the estimator of gjo is 

&(*) = 22 W 0im,h j {t){Y jim - X$ im ij - (gj(tjim) ~ 9j(t))}- 
i=l m=l 

Based on this modified estimator <jj(t), we define the auxiliary variable 



T 

±Ti KnimiGj) V h 3 J 



=1 "Jtmydj) 

x {Y jim - X T jim ij - gj (t) - (g~j(t jim ) - 9j(t))} 
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for empirical likelihood formulation. At true function gjo(t), E[Rji{gjo(t)}] = 
o(l). 

Using a similar procedure to L nj ((3j) as given in (3.5) and (3.6), the 
empirical likelihood for gjo(t) is 

L n 3 {9jo(t)} = maxj JJpjij 

subject to Y^LiPji = 1 an d Y^LiPjiRji{9j{^)} = 0- The latter is obtained 
in a similar fashion as we obtain (3.5) by introducing Lagrange multipliers 
so that 



M&o(i)}=n{-^ r 

i=l ^ 3 



+ i] j (t)R ji {g j0 (t)} J ' 
where r]j{t) is a Lagrange multiplier that satisfies 

(4.3) T R ^ o{t)} = 0. 

The log empirical likelihood ratio for gio{t) = • • • = gko(t) ■= g(t), say, is 

k nj 

(4.4) C n {t) = 2 min ^ ^ fog( 1 + m (t)Rji{g(t)}), 

3 =i i=i 

which is analogues of t n in (3.7). As shown in the proof of Theorem 3 given 
in the supplement article [Chen and Zhong (2010)], the leading order term 
of the C n {t) is a studentized version of the distance 

(9i(t) - 92(t), 9l (t) -&(*),.- -,9i(t) - g k (t)), 

namely between g±(t) and the other gj(t)(j ^ 1). This motivates us to pro- 
pose using 

(4.5) T n = [ C n {t)w{t)dt 

Jo 

to test for the equivalence of {gjo(")}j=ii where w(t) is a probability weight 
function over [0, 1]. 

To define the asymptotic distribution of 7~ n , we assume without loss of 
generality that for each hj and Tj, j = 1, . . . , k, there exist fixed finite positive 
constants a.j and bj such that ctjTj = T and bjhj = h for some T and h as 
h — > 0. Effectively, T is the smallest common multiple of Ti,...,Tfc. Let 
K^(t) = fK{w)K(t - cw)dt and ^ 4) (0) = f K^ 2 \wy/c)K[f c (w/y/c) dw. 
For c= 1, we resort to the standard notations of K^- 2 \t) and ET( 4 )(0) for 
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K^\t) and K^\o), respectively. For each treatment j, let fj be the super- 
population density of the design points {tji m }. Let aj = p~ x a.j, 



and V 3 {t) = K^{Q)ay 3 {t) where a% = ^ E2=i E i^^)}- 
thermore, we define 

A(i) = ^67 1 /^ (4) (0)(l-^(t)) 2 



and 



Ml 



-3=1 

We consider a sequence of local alternative hypotheses: 

(4.6) g j0 (t) = gio (t) + C jn A nj (t), 

where C in = (njTj)- 1 / 2 ^ 1 ^ for j = 2, . . . ,k and {A nj (t)} n >i is a sequence 
of uniformly bounded functions. 

Theorem 3. Assume conditions A1-A4 in the Appendix and h = O^" 1 ^), 
then under (4-6), 

h- l / 2 (T n -Ho)^N(0,o- 2 ), 
where fi = (k - 1) + h^ 2 m and a 2 = 2A'( 2 )(0)" 2 A(t)m 2 (t) dt. 

We note that under i? 6 : 9io(') = ■ ■• = <7fco(-)> ^nj(t) = which yields 
Hi = and 

h- l ' 2 {T n - (k -l)}^N(0,a 2 ). 

This may lead to an asymptotic test at a nominal significance level a that 
rejects if 

(4.7) T n >h 1 / 2 a z a + (k-l), 



w(t) dt. 
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where z a is the upper a quantile of iV(0, 1) and cto is a consistent estimator 
of <jq. The asymptotic power of the test under the local alternatives is 1 — 
&(z a — ^-), where $(•) is the standard normal distribution function. This 
indicates that the test is powerful in differentiating null hypothesis and its 
local alternative at the convergence rate 0{n- for Cj n . The rate is 

the best when a single bandwidth is used [Hardle and Mammen (1993)]. 

If all the hj(j = 1, . . . , k) are the same, the asymptotic variance Cq = 2(k — 
1)K& {0)- 2 K^(0) Jo 1 w 2 (t) dt, which means that the test statistic under 
Hq}, is asymptotic pivotal. However, when the bandwidths are not the same, 
which is most likely as different treatments may require different amount 
of smoothness in the estimation of gjo(-), the asymptotical pivotalness of 
T n is no longer available, and estimation of a 1 is needed for conducting 
the asymptotic test in (4.7). We will propose a test based on a bootstrap 
calibration to the distribution of T n in Section 6. 

Remark 4. Similar to Remarks 1 and 2 made on the ANOVA tests 
for the covariate effects, the proposed ANOVA test for the nonpar ametric 
baseline functions (Theorem 3) remains valid in the absence of missing values 
or if the missing propensity is misspecified as long as the responses do not 
contribute to the missingness. 

Remark 5. We note that the proposed test is not affected by the within- 
subject dependent structure (the longitudinal aspect) due to the fact that 
the formulation of the empirical likelihood is made for each subject. This 
is clearly shown in the construction of Rji{gj(t)} and by the fact that the 
nonparametric functions can be separated from the covariate effects in the 
semiparametric model. Again this would be changed if we are interested in 
estimation as the correlation structure in the longitudinal data will affect 
the estimation efficiency. However, the test will be dependent on the choice 
of the weight function w(-), and {aj}, {pj} and {bj}, the relative ratios 
among {!}}, {nf\ and {hj}. 

Remark 6 . The ANOVA test statistics for the time effects for the semi- 
parametric model can be readily modified to obtain ANOVA test for purely 
nonparametric regression by simply setting £j = in the formulation of the 
test statistic C n (t). In this case, the model (2.1) takes the form 

(4.8) Yji{t) = gj (Xji(t),t) + eji(t), 

where gj(-) is the unknown nonparametric function of Xji(t) and t. The 
proposed ANOVA test can be viewed as generalization of the tests considered 
in Mund and Dettle (1998), Pardo- Fernandez, Van Keilegom and Gonzalez- 
Manteiga (2007) and Wang, Akritas and Van Keilegom (2008) by considering 
both the longitudinal and missing aspects. See also Cao and Van Keilegom 
(2006) for a two sample test for the equivalence of two probability densities. 
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5. Tests on interactions. Model (1.1) contains an interactive term 
M(Xji m ,t) that is flexible in prescribing the interact between Xji m and 
the time, as long as the positive definite condition in condition A3 is sat- 
isfied. In this section, we propose tests for the presence of the interaction 
in the jth treatment and the ANOVA hypothesis on the equivalence of the 
interactions among the treatments. 

We firstly consider testing Hq c : 7^0 = vs. H\ c : 7^0 7^ for a fixed j. In the 
formulation of the empirical likelihood for 7^0, we treat Mji m = M(Xji m ,t) 
as a covariates with the same role like Xji m in the previous section when 
we constructed empirical likelihood for /3jQ. For this purpose, we define es- 
timating equations for 7^0 

(5.1) 0ji(7jo) = Yl — 3J 7Z^M jim {Y jim - XJ im (3j - MJ im ~/ j0 ), 



where 



(5.2) 



[ 5- 



,i=l m=l Kjimyvj, 



^ ^ Jl ™z Xjim(Yjim - ^jimljQ} 



i=l m=l n 3im 



9j) 



is the "estimator" of (3j at the true 7^0- Similar to establishing l n j(Pj) in 
Section 3, the log-empirical likelihood for 7,0 can be written as 



^(7i) = 2£log{l + A^ji(7i)}, 

i=l 



where the Lagrange multipliers Aj satisfies 

tain) 



M El 



-1 + A^( 7i ) 



0. 



To test for H od : 7l0 = 720 = • • • = 7 fc0 vs. H ld : 7 i0 / 7jo for some i^j, we 
construct the joint empirical likelihood ratio 

k rij 

(5.4) P n := 2minEE 1 °g{ 1 + A J<Mt)}, 

7 J'=l i=l 

where Aj satisfy (5.3). 

The asymptotic distributions of the empirical likelihood ratios . (0) and 
In under the null hypotheses are given in the next theorem whose proofs 
will not be given as they follow the same routes in the proof of Theorem 1 
by exchanging Xji m and /3jo with Mji m and jjo, respectively. 
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Theorem 4. Under conditions A1-A4 given in the Appendix, then (i) un- 
der H 0c , llj (0) A \ 2 q as rij -> oo; (ii) under H od , ll -4- xL-l)q as n ~^ 00 • 

Based on Theorem 4, an a-level empirical likelihood ratio test for the pres- 
ence of the interaction in the jth sample rejects Hq c if £nj(0) > x\ a i anc ^ 
the ANOVA test for the equivalence of the interactive effects rejects Hq^ if 
Gi > X{j~-i) qa - The ANOVA test for has a similar local power perfor- 
mance as that described after Theorem 2 for the ANOVA test regarding /3jo 
in Section 3. The power properties of the test for H$ c can be established 
using a much easier method. 

We have assumed parametric models for the interaction in model (1.1). 
A semiparametric model would be employed to model the interaction given 
that the model for the time effect is nonparametric. The parametric in- 
teraction is a simplification and avoids some of the involved technicalities 
associated with a semiparametric model. 

6. Bootstrap calibration. To avoid direct estimation of Cg in Theorem 3 
and to speed up the convergence of T n , we resort to the bootstrap. While 
the wild bootstrap [Wu (1986), Liu (1988) and Hardle and Mammen (1993)] 
originally proposed for parametric regression without missing values has 
been modified by Shao and Sitter (1996) to take into account missing values, 
we extend it further to suit the longitudinal feature. 

Let tj and vj 1 be the sets of the time points with full and missing obser- 
vations, respectively. According to model (2.2), we impute a missing Xji(t) 
from Xji(t),t £ t°, so that for any tEvj 1 

(6.1) Xji(t) = ^2 Wji m , hj (t)X jim , 

i=l m=l 

where 10^^.(4) is the kernel weight defined in (3.2). 

To mimic the heteroscedastic and correlation structure in the longitudinal 
data, we estimate the covariance matrix for each subject in each treatment. 
Let 

An estimator of cj(i), the variance of £ji(t), is o^{t) = J2i=i Em=i w jim,hj(t) x 
ij im and an estimator of pj(s,t), the correlation coefficient between £ji(t) 
and Sji(s) for s / 1, is 

Pj(s, i) = y y Hji m m i(s,t^eji m eji m i, 
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where Cjirrt — £jim/ 1 @ jitjim) j 

y ^ ^ 3jirn3jim'K-i,j(s tjim)^bj\P tjim , )/'Kjim,m'(@j^) 

Sj=l Em^m' fijimfijim' K-bj — tjim^jK^. (t tjim')/'Kjim,m'{@j) 

and iTji m>m >{0j) = TTjim{0j)TT jim >(9j) if |m - m'| > d; ^ji m , m '{Qj) = ^jim b (0j) 
\{ \m — m!\ < d where = max(m,m'). Here 6j is a smoothing band- 
width which may be different from the bandwidth hj for calculating the 
test statistics T n [Fan, Huang and Li (2007)]. Then, the covariance Sjj of 
£ji = (£jn, ■ • ■ , £jiTj) T is estimated by Sjj which has aj(tji m ) as its mth diag- 
onal element and Pj{tjik,tjii)aj{tjik)(jj{tjii) as its (fc, Z)th element for k^l. 

Let Yji,5ji,tji be the vector of random variables of the (j, i)th subject, 
Xji = (Xji(tjn), . . . ,Xj i (tjiT j )) r and gjo{t s i) = {gjo{t s n), ... ,gjo{t s iT k )) T , 
where s may be different from j. Let = where X^ contains 

observed Xji(t) for tj G t° and X™ collects the imputed Xji(t) for tef™ 
according to (6.1). Plugging the value of X^, we get Mj^ = {Mj^Mj^ 1 }, the 
observed and the imputed interactions for (j,z)th subject and then X^. 
The proposed bootstrap procedure consists of the following steps: 
Step 1. Generate a bootstrap re-sample {Y^X^,^,^} for the (j,i)th 
subject by 

Yji = ^ji T ij + 9i {tji) + £jie*ji, 
where e*j's are i.i.d. random vectors simulated from a distribution satisfying 

E{e* i ) = and Var(e*j) = J^, S* im ~ Bernoulli(7Tjj m (^)) where Oj is esti- 
mated based on the original sample as given in (2.5). Here, g\{tji) is used as 
the common nonparametric time effect to mimic the null hypothesis Hq^. 

Step 2. For each treatment j, we reestimate £j, 8j and gj(t) based on 
the resample {Y* i ,'X < j i ,5* i ,tji\ and denote them as £*, 6* and g*j{t). The 
bootstrap version of Rji{gi(t)} is 

x {Y* im - X^J* - gi (t) - {g*(t jm ) - g*(t)}} 

and use it to substitute Rji{gj(t)} in the formulation of C n (t), we obtain 
£*(*) and then T* = J C* n {t)w{t) dt. 

Step 3. Repeat the above two steps B times for a large integer B and 
obtain B bootstrap values {T* b } b=l . Let i a be the 1 — a quantile of {T* b }^ =l , 
which is a bootstrap estimate of the 1 — a quantile of T n - Then, we reject 
the null hypothesis Hqi, if T n > t a . 

The following theorem justifies the bootstrap procedure. 
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Theorem 5. Assume conditions A1-A4 in the Appendix hold and h = 
0(?i -1 / 5 ). Let X n denote the original sample, h and a\ he defined as in The- 
orem 3. The conditional distribution of h~ x l 2 {J~* — uq) given X n converges 
to iV(0,(To) almost surely, namely, 

h~ 1/2 {V -(k- l)}\X n A N(0, (j 2 ) a.8. 

7. Simulation results. In this section, we report results from simulation 
studies which were designed to confirm the proposed ANOVA tests pro- 
posed in the previous sections. We simulated data from the following three- 
treatment model: 

Yji m — Xji m f3j + Mji rn '~ % /j -\- 9j(tjirn) ~i~ &jim and 

(7.1) 

Xji m — 2 1.5tjfj m -|- Uji m , 

where M jim = t j i rn x (X jim - 1.5) 2 , e jim = eji + u jim , u jim ~ N(0,al.), eji ~ 
N(0, cr 2 .) and Uj im ~ N(0, c 2 .) for j = {1, 2, 3}, i = 1, ... ,nj and m = 1, . . . ,Tj. 

This structure used to generate {eji m } r ^ =1 ensures dependence among the 
repeated measurements {Yji m } for each subject i. The correlation between 

Yjim and Yju for any m ^ I is o% /(of + cr 2 .). The time points {i,j m }^ =1 
were obtained by first independently generating uniform[0, 1] random vari- 
ables and then sorted in the ascending order. We set the number of repeated 
measures Tj to be the same, say T, for all three treatments; and chose T = 5 
and 10, respectively. The standard deviation parameters in (7.1) were a ai = 
0.5,0"^ = 0.5, a Cl = 0.2 for the first treatment, a a2 = 0.5,<7ft 2 = 0.5, a C2 = 0.2 
for the second and a a3 =0.6, o"b 3 =0.6, a Cz =0.3 for the third. 

The parameters and the time effects for the three treatments were: 



Treatment 1 
Treatment 2 
Treatment 3 



/3 1 =2, 7l = l, 5l (t)=2sin(27rt); 

fa = 2 + £> 2n ,72 = 1 + D 2n ,g 2 (t) = 2sin(2vrt) - A 2 „(t); 
/? 3 = 2 + D 3n ,73 = 1 + D 3n ,g 3 (t) = 2sin(27rt) - A 3n (i). 



We designated different values of £>2n> i^3n, ^2n(t) and A3„(i) in the evalu- 
ation of the size and the power, whose details will be reported shortly. 

We considered two missing data mechanisms. In the first mechanism (I), 
the missing propensity was 

(7.2) logit{P(<5j im = l|(5j imim _i = l } Xji,Yji)} = 9jXj^ m ^ for m > 1, 

which is not dependent on the response Y, with 9\ = 3,02 = 2 and #3 = 2. 
In the second mechanism (II), 

logit{P(5 jim = l|<5j- imm _i = l,Xji,Yji)} 

( 7 - 3 ) 

9jl X ji(m-l) + 6j2{Yji(m-l) - Y ji(m-2)}, if m > 2, 

^ji^ii( m _i), ifm = 2; 



18 



S. X. CHEN AND P.-S. ZHONG 



which is influenced by both covariate and response, with 6\ = (0u,6i 2 ) T = 
(2, -l) r ,6 2 = (6 2l ,6 22 y = (2,-1.5) T and 6 3 = (#31, # 32 ) T = (2, -1-5) T . In both 
mechanisms, the first observation (m = 1) for each subject was always ob- 
served as we have assumed earlier. 

We used the Epanechnikov kernel K(u) = 0.75(1 — ii 2 )+ throughout the 
simulation where (•)+ stands for the positive part of a function. The band- 
widths were chosen by the "leave-one-subject" out cross-validation. Specifi- 
cally, we chose the bandwidth hj that minimized the cross-validation score 
functions 

\ V °i im (Y -X T - B { ~ i) -M T - - a { ~ i] (t-- )) 2 

i=l m=l Kjimivj) 

where , 7: l ' and g~- (tji m ) were the corresponding estimates with- 
out using observations of the ith. subject. The cross-validation was used to 
choose an optimal bandwidth for representative data sets and fixed the cho- 
sen bandwidths in the simulations with the same sample size. We fixed the 
number of simulations to be 500. 

The average missing percentages based on 500 simulations for the missing 
mechanism I were 8%, 15% and 17% for treatments 1-3, respectively, when 
T = 5, and were 16%, 28% and 31% when T = 10. In the missing mecha- 
nism II, the average missing percentages were 10%, 8% and 15% for T = 5, 
and 23%, 20% and 36% for T = 10, respectively. 

For the ANOVA test for H$ a : /3\q = (3 2 q = /330 with respect to the covariate 
effects, three values of D 2n and D^ n : 0, 0.2 and 0.3, were used, respectively, 
while A 2n (t) = A^ n (t) = 0. Table 1 summarizes the empirical size and power 
of the proposed EL ANOVA test with 5% nominal significant level for Ho a 
for 9 combinations of (D 2n , D^ n ), where the sizes corresponding to D 2n = 
and D^ n = 0. We observed that the size of the ANOVA tests improved as the 
sample sizes and the observational length T increased, and the overall level 
of size were close to the nominal 5%. This is quite reassuring considering 
the ANOVA test is based on the asymptotic chi-square distribution. We 
also observed that the power of the test increased as sample sizes and T 
were increased, and as the distance among the three /3jo was increased. For 
example, when D 2n = 0.0 and D% n = 0.3, the L 2 distance was V0.3 2 + 0.3 2 = 
0.424, which is larger than VO.l 2 + 0.2 2 +0.3 2 = 0.374 for D 2n = 0.2 and 
Z?3„ = 0.3. This explains why the ANOVA test was more powerful for D 2n = 
0.0 and D% n = 0.3 than D 2n = 0.2 and D^ n = 0.3. At the same time, we see 
similar power performance between the two missing mechanisms. 

To gain information on the empirical performance of the test on the ex- 
istence of interaction, we carried out a test for Hq c : 720 = 0. In the sim- 
ulation, we chose 720 = 0,0.2,0.3,0.4, /3 2 q = 2 + 720 and fixed A 2n (t) = 0, 
respectively. Table 2 summarizes the sizes and the powers of the test. Ta- 
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Table 1 



Empirical size 


and 


power of the 5% ANOVA test f, 


nr Ho a 


:/3io = 


P20 = 030 




Sample size 


D- 2n 




Missingness 


T 


Missingness 


n± n 2 n 3 


T I 


II 


I 


II 


60 65 55 


0.0 


0.0 (size) 


5 0.042 


0.050 


10 


0.046 


0.044 




0.2 


0.0 


0.192 


0.254 




0.408 


0.434 




0.3 


0.0 


0.548 


0.630 




0.810 


0.864 




n n 

u.u 


0.2 


0.236 


0.214 




0.344 


0.354 




0.0 


0.3 


0.508 


0.546 




0.714 


0.722 




0.2 


0.2 


0.208 


0.262 




0.446 


0.458 




0.2 


0.3 


0.412 


0.440 




0.680 


0.698 




0.3 


0.2 


0.426 


0.490 




0.728 


0.728 




0.3 


0.3 


0.594 


0.620 




0.836 


0.818 


100 110 105 


0.0 


0.0 (size) 


5 0.052 


0.054 


10 


0.042 


0.038 




0.2 


0.0 


0.426 


0.470 




0.686 


0.718 




0.3 


0.0 


0.854 


0.854 




0.964 


0.974 




0.0 


0.2 


0.406 


0.444 




0.612 


0.568 




0.0 


0.3 


0.816 


0.836 




0.936 


0.910 




0.2 


0.2 


0.404 


0.480 




0.674 


0.686 




0.2 


0.3 


0.744 


0.694 




0.944 


0.882 




0.3 


0.2 


0.712 


0.768 




0.922 


0.920 




0.3 


0.3 


0.824 


0.814 




0.972 


0.970 



ble 3 reports the simulation results of the ANOVA test on the interaction 
effect Hod : 710 = 720 = 730 with a similar configurations as those used as the 
ANOVA tests for the covarites effects reported in Table 1. We observe satis- 
factory performance of these two tests in terms of both the accurate of the 
size approximation and the empirical power. In particular, the performance 
of the ANOVA tests were very much similar to that conveyed in Table 1. 

Table 2 

Empirical size and power of the 5% test for the existence of interaction Hoc ■ 720 = 



Sample size Missingness Missingness 



m 


n 2 


n 3 


720 


T 


I 


II 


T 


I 


II 


60 


65 


55 


0.0 (size) 


5 


0.052 


0.048 


10 


0.048 


0.052 








0.2 




0.428 


0.456 




0.568 


0.636 








0.3 




0.722 


0.788 




0.848 


0.882 








0.4 




0.928 


0.952 




0.948 


0.968 


100 


110 


105 


0.0 (size) 


5 


0.054 


0.046 


10 


0.056 


0.042 








0.2 




0.608 


0.718 




0.694 


0.812 








0.3 




0.940 


0.938 




0.940 


0.958 








0.4 




0.986 


0.994 




0.952 


0.966 
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Table 3 

Empirical size and power of the 5% ANOVA test for H(,d :71c = 720 = 730 



Sample size Missingness Missingness 



m 


»2 


n 3 


D 2n 


D 3n 


T 


I 


II 


T 


I 


II 


60 


65 


55 


0.0 


0.0 (size) 


5 


0.058 


0.058 


10 


0.068 


0.036 








0.2 


0.0 




0.134 


0.188 




0.232 


0.254 








0.3 


0.0 




0.358 


0.486 




0.510 


0.622 








0.0 


0.2 




0.136 


0.166 




0.230 


0.218 








0.0 


0.3 




0.356 


0.414 




0.466 


0.474 








0.2 


0.2 




0.170 


0.208 




0.286 


0.276 








0.2 


0.3 




0.292 


0.328 




0.462 


0.428 








0.3 


0.2 




0.266 


0.356 




0.498 


0.474 








0.3 


0.3 




0.392 


0.476 




0.578 


0.588 


100 


110 


105 


0.0 


0.0 (size) 


5 


0.068 


0.040 


10 


0.054 


0.046 








0.2 


0.0 




0.262 


0.366 




0.354 


0.432 








0.3 


0.0 




0.654 


0.744 




0.744 


0.820 








0.0 


0.2 




0.272 


0.330 




0.340 


0.334 








0.0 


0.3 




0.590 


0.676 




0.722 


0.672 








0.2 


0.2 




0.282 


0.332 




0.412 


0.410 








0.2 


0.3 




0.528 


0.582 




0.716 


0.640 








0.3 


0.2 




0.502 


0.580 




0.680 


0.728 








0.3 


0.3 




0.672 


0.674 




0.814 


0.808 



We then evaluate the power and size of the proposed ANOVA test regard- 
ing the nonparametric components. To study the local power of the test, 
we set /S.2n{t) = U n sin(27rt) and A^ n (t) = 2sin(2-7rf) — 2sin(27r(t + V n )), and 
fixed D2n = and Ds n = 0.2. Here, U n and V n were designed to adjust the 
amplitude and phase of the sine function. The same kernel and bandwidths 
chosen by the cross-validation as outlined earlier in the parametric ANOVA 
test were used in the test for the nonparametric time effects. We calculated 
the test statistic T n with w{t) being the kernel density estimate based on all 
the time points in all treatments. We applied the wild bootstrap proposed 
in Section 6 with B = 100 to obtain to. 05) the bootstrap estimator of the 5% 
critical value. The simulation results of the nonparametric ANOVA test for 
the time effects are given in Table 4. 

The sizes of the nonparametric ANOVA test were obtained when U n = 
and V n = 0, which were quite close to the nominal 5%. We observe that the 
power of the test increased when the distance among <7i(-), gi{-) and g${-) 
were becoming larger, and when the sample size or repeated measurement 
T were increased. We noticed that the power was more sensitive to change 
in V n , the initial phase of the sine function, than U n . 

We then compared the proposed tests with a test proposed by Scheike 
and Zhang (1998). Scheike and Zhang's test was comparing two treatments 
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Table 4 

Empirical size and power of the 5% ANOVA test for : gi(-) = <?2(-) = fl3(-) with 
A 2 „(t) = U„sm(2TTt) and A 3n (t) = 2sin(27rf) - 2sin(27r(t + V n )) 



Sample size Missingness Missingness 



m 


n-2 


n 3 


u n 


V n 


T 


I 


II 


T 


I 


II 


60 


65 


55 


0.00 


0.00 (size) 


5 


0.040 


0.050 


10 


0.054 


0.060 








0.30 


0.00 




0.186 


0.232 




0.282 


0.256 








0.50 


0.00 




0.666 


0.718 




0.828 


0.840 








0.00 


0.05 




0.664 


0.726 




0.848 


0.842 








0.00 


0.10 




1.000 


1.000 




1.000 


1.000 


100 


110 


105 


0.00 


0.00 (size) 


5 


0.032 


0.062 


10 


0.050 


0.036 








0.30 


0.00 




0.434 


0.518 




0.526 


0.540 








0.50 


0.00 




0.938 


0.980 




0.992 


0.998 








0.00 


0.05 




0.916 


0.974 




1.000 


1.000 








0.00 


0.10 




1.000 


1.000 




1.000 


1.000 



for the nonparametric regression model (4.8) for longitudinal data without 
missing values. Their test was based on a cumulative statistic 

T(z)= [ (h{t)-fo{t))dt, 

where a, z are in a common time interval [0, 1]. They showed that \Jn\ + n 2 T(z) 
converges to a Gaussian Martingale with mean and variance function 
p^ 1 h\(z) + p^ 1 h 2 {z) , wherehj(z) = (Tj(y) fj l (y) dy . Hence, the test statis- 
tic T(l — a)/ y / 'var{T(l — a)} is used for two group time-effect functions 
comparison. 

To make the proposed test and the test of Scheike and Zhang (1998) 
comparable, we conducted simulation in a set-up that mimics the setting 
of model (7.1) but with only the first two treatments, no missing values 
and only the nonparametric part in the regression by setting /?,- = 7j =0. 
Specifically, we test for H :g 1 (-) = g 2 {-) vs. H 1 :g 1 (-) = g 2 (-) + A 2n (-) for 
three cases of the alternative shift function A2 n ( - ) functions which are spelt 
out in Table 5 and set a = in the test of Scheike and Zhang. The simulation 
results are summarized in Table 5. We found that in the first two cases (I 
and II) of the alternative shift function A2 n , the test of Scheike and Zhang 
had little power. It was only in the third case (III), the test started to pick 
up some power although it was still not as powerful as the proposed test. 

8. Analysis on HIV-CD4 data. In this section, we analyzed a longitu- 
dinal data set from AIDS Clinical Trial Group 193A Study [Henry et al. 
(1998)], which was a randomized, double-blind study of HIV- AIDS patients 
with advanced immune suppression. The study was carried out in 1993 with 
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Table 5 

The empirical sizes and powers of the proposed test (CZ) and the test (SZ) proposed by 
Scheike and Zhang (1998) for Hob :<?i(-) = ff2(-) vs. Hn -gi(-) = <?2(-) + A2n(-) 



Sample size Tests Tests 



Til 


Tl2 


W'3 


77 T 


CZ sz 


T 


CZ 


SZ 


60 


65 


55 




Case I: A 2n {t) = U n 


sin(27rf) 












0.00 (size) 5 


0.060 0.032 


10 


0.056 


0.028 








0.30 


0.736 0.046 




0.844 


0.028 








u.ou 


1 nnn n n/is 

l.UUU U.U4o 




1 nnn 


U.UZD 








Case II: 


A 2n (t) = 2sin(27rf) - 


2sin(27r 


(t + V n )) 










0.05 


1.000 0.026 




1.000 


0.042 








0.10 


1.000 0.024 
Case III: A 2 „(t) = 


= -U n 


1.000 


0.044 








0.10 


0.196 0.162 




0.206 


0.144 








0.20 


0.562 0.514 




0.616 


0.532 


100 


110 


105 




Case I: A2„(t) = U n 


sin(27rt) 












0.00 (size) 5 


0.056 0.028 


10 


0.042 


0.018 








0.30 


0.982 0.038 




0.994 


0.040 








0.50 


1.000 0.054 




1.000 


0.028 








Case II: 


A 2n (t) = 2sin(27rf) - 


2sin(27r 


(t+u n )) 










0.05 


1.000 0.022 




1.000 


0.030 








0.10 


1.000 0.026 
Case III: A 2n (t) = 


= -u n 


1.000 


0.030 








0.10 


0.290 0.260 




0.294 


0.218 








0.20 


0.780 0.774 




0.760 


0.730 



1309 patients who were randomized to four treatments with regard to HIV- 
1 reverse transcriptase inhibitors. Patients were randomly assigned to one 
of four daily treatment regimes: 600 mg of zidovudine alternating monthly 
with 400 mg didanosine (treatment I); 600 mg of zidovudine plus 2.25 mg of 
zalcitabine (treatment II); 600 mg of zidovudine plus 400 mg of didanosine 
(treatment III) ; or 600 mg of zidovudine plus 400 mg of didanosine plus 400 
mg of nevirapine (treatment VI). The four treatments had 325, 324, 330 and 
330 patients, respectively. 

The aim of our analysis was to compare the effects of age (Age) , baseline 
CD4 counts (PreCD4) and gender (Gender) on 7= log(CD4 counts +1). 
The semiparametric model regression is, for j = 1,2,3 and 4, 

(8.1) Yji (t) = p n Ag eji +f3 j2 PreCD4 ji +/3 i3 Gender^ + 9j (t) + Eji (t) 

with the intercepts absorbed in the nonparametric gj(-) functions, and j3j = 
(fiji, /3j2, fijzY is the regression coefficients to the covariates (Age, PreCD4, 
Gender). 
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To make gj(t) more interpretable, we centralized Age and PreCD4 so that 
their sample means in each treatment were 0, respectively. As a result, gj(t) 
can be interpreted as the baseline evolution of Y for a female (Gender = 0) 
with the average PreCD4 counts and the average age in treatment j. This 
kind of normalization is used in Wu and Chiang (2000) in their analyzes 
for another CD4 data set. Our objectives were to detect any difference in 
the treatments with respect to (i) the covariates; and (ii) the nonparametric 
baseline functions. 

Measurements of CD4 counts were scheduled at the start time 1 and at 
a 8-week intervals during the follow-up. However, the data were unbalanced 
due to variations from the planned measurement time and missing values 
resulted from skipped visits and dropouts. The number of CD4 measure- 
ments for patients during the first 40 weeks of follow-up varied from 1 to 
9, with a median of 4. There were 5036 complete measurements of CD4, 
and 2826 scheduled measurements were missing. Hence, considering missing 
values is very important in this analysis. Most of the missing values follow 
the monotone pattern. Therefore, we model the missing mechanism under 
the monotone assumption. 

We considered three logistic regression models for the missing propensities 
and used the AIC and BIC criteria to select the one that was the mostly 
supported by data. The first model (Ml) was a logistic regression model 
for Pj(Xjit t 3,Yjit t 3\ 9jo) that effectively depends on Xy lt (the PreCD4) and 
(Yjj( t _i), Yji( t -2),Yji(t-3)) if t > 3. For t < 3, it relies on all Y jit observed 
before t. In the second model (M2), we replace the Xjn in the first model 
with an intercept. In the third model (M3), we added to the second logistic 
model with covariates representing the square of Y^( t _i) and the interactions 
between Yju^\\ and Yj i ( t _ 2 y In the formulation of the AlC and BIC criteria, 
we used the binary conditional likelihood given in (2.5) with the respective 
penalties. The difference of AIC and BIC values among these models for four 
treatment groups is given in Table 6. Under the BIC criterion, M2 was the 
best model for all four treatments. For treatments II and III, M3 had smaller 
AIC values than M2, but the differences were very small. For treatments I 
and VI, M2 had smaller AIC than M3. As the AIC tends to select more 
explanatory variables, we chose M2 as the model for the parametric missing 
propensity. 

Model (8.1) does not have interactions. It is interesting to check if there 
is an interaction between gender and time. Then the model becomes 

Yji(t) = + (3j2 PreCD4 ii + p j3 Gender^ 

(8.2) 

+ 7j4 Gender x t + gj (t) + Eji{t). 

We applied the proposed test in Section 5 for HQ C :^j^ = for j = 1,2,3 
and 4, respectively. The p- values were 0.9234,0.9885,0.9862 and 0.5558, re- 
spectively, which means that the interaction was not significant. Therefore, 
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Table 6 

Difference in the AIC and BIC scores among three models (}A1)-(M3) 





Treatment I 


Treatment II 


Treatment III 


Treatment VI 


Models 


AIC BIC 


AIC BIC 


AIC BIC 


AIC BIC 


(M1)-(M2) 
(M2)-(M3) 


3.85 3.85 
-2.47 -11.47 


14.90 14.90 
0.93 -8.12 


17.91 17.91 
0.30 -8.75 


10.35 10.35 
-3.15 -12.27 



in the following analyzes, we would not include the interaction term and 
continue to use model (8.1). 

Table 7 reports the parameter estimates f3j of f3j based on the estimating 
function Zji(fij) given in Section 3. It contains the standard errors of the 
estimates, which were obtained from the length of the EL confidence inter- 
vals based on the marginal empirical likelihood ratio for each j3j as proposed 
in Chen and Hall (1994). In getting these estimates, we use the "leave-one- 
subject" cross-validation [Rice and Silverman (1991)] to select the smoothing 
bandwidths {hj}j =l for the four treatments, which were 12.90, 7.61, 8.27 and 
16.20, respectively. We see that the estimates of the coefficients for the Age 
and PreCD4 were similar among all four treatments with comparable stan- 
dard errors, respectively. In particular, the estimates of the Age coefficients 
endured large variations while the estimates of the PreCD4 coefficients were 
quite accurate. However, estimates of the Gender coefficients had different 
signs among the treatments. We may also notice that the confidence intervals 
from treatments I— IV for each coefficient were overlap. 

We then formally tested Ho a : f}\ = 02 = 03 = 04- The empirical likelihood 
ratio statistic l n was 8.1348, which was smaller than Xg o 95 = 16-9190, which 
produced a p-value of 0.5206. So we do not have enough evidence to reject 
Ho a at a significant level 5 %. The parameter estimates reported in Table 7 
suggested similar covariate effects between treatments I and II, and between 
treatments III and IV, respectively; but different effects between the first two 
treatments and the last two treatments. To verify this suggestion, we carry 
out formal ANOVA test for pair-wise equality among the 0j's as well as for 
equality of any three /3j's. The p- values of these ANOVA test are reported in 



Table 7 

Parameter estimates and their standard errors 





Treatment I 


Treatment II 


Treatment III 


Treatment IV 


Coefficients 


01 


02 


03 


04 


Age 

PreCD4 
Gender 


0.0063 (0.0039) 
0.7308 (0.0462) 
0.1009 (0.0925) 


0.0050 (0.0040) 
0.7724 (0.0378) 
0.1045 (0.0920) 


0.0047 (0.0058) 
0.7587 (0.0523) 
-0.3300 (0.1510) 


0.0056 (0.0046) 
0.8431 (0.0425) 
-0.3055 (0.1136) 
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Table 8 
p-values of ANOVA tests for fij 's 



Hoa 


p- value 


Ho a 


p- value 


Px=p2 


0.9661 


Pl=P2= Ps 


0.7399 


Pi = Ps 


0.4488 


Pl=P 2 = Pi 


0.4011 


Pl=Pi 


0.1642 


Pi=P 3 = Pi 


0.3846 


P 2 =Ps 


0.4332 


P2=PA=P4 


0.4904 


132=134 


0.2523 


Pl=P2=P3= Pi 


0.5206 


P 3 =Pi 


0.8450 







Table 8. Indeed, the difference between the first two treatments and between 
the last two treatments were insignificant. However, the differences between 
the first three (I, II and III) treatments and the last treatment were also not 
significant. 

We then tested for the nonpar ametric baseline time effects. The kernel 
estimates gjit) are displayed in Figure 1, which shows that treatments I and 
II and treatments III and IV had similar baselines evolution overtime, re- 
spectively. However, a big difference existed between the first two treatments 
and the last two treatments. Treatment IV decreased more slowly than that 
of the other three treatments, which seemed to be the most effective in slow- 
ing down the decline of CD4. We also found that during the first 16 weeks 
the CD4 counts decrease slowly and then the decline became faster after 16 
weeks for treatments I, II and III. 



Treatment I Treatment II 




Week(t) Week(t) Week(t) 



(a) (b) 

Fig. 1. (a) The raw data excluding missing values plots with the estimates of gj(t) 
(j = 1,2,3,4:). (b) The estimates of gj(t) in the same plot: treatment I (solid line), treat- 
ment II (short dashed line), treatment III (dashed and doted line) and treatment IV (long 
dashed line). 
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Table 9 

p-values of A NOVA tests on gj(-)'s 



Hob 


p- value 


Hob 


p- value 


Sl(') = 32(0 


0.894 


31 (•) =32(0 =33(0 


0.046 


ffi(0=5s(0 


0.018 


ffl(-) =32(0 =34(0 


0.010 


ffi(0 =34(0 


0.004 


3i(0 =33(0 =34(0 


0.000 


32(-) =3a(-) 


0.020 


32(0 = 33(0 = 34(0 


0.014 


32(0 = 34(0 


0.006 


3i(0 =32(0 =33(0 =34(0 


0.004 


33(0 = 34(0 


0.860 







The p-value for testing Ho& : 5i(") = 52( - ) = 53Q = 54(0 is shown in Ta- 
ble 9. The entries were based on 500 bootstrapped resamples according to the 
procedure introduced in Section 6. The statistics T n for testing ifofc : 5i(") = 
52(0 = 53(") = 54(") was 3965.00, where we take w(t) = 1 over the range of t. 
The p- value of the test was 0.004. Thus, there existed significant difference in 
the baseline time effects 5j(-)' s among treatments I-IV. At the same time, we 
also calculate the test statistics T n for testing gi(-) = 52(0 and gs(-) = #4(-)- 
The statistics values were 19.26 and 26.22, with p-values 0.894 and 0.860, 
respectively. These p-values are much bigger than 0.05. We conclude that 
treatment I and II has similar baseline time effects, but they are significantly 
distinct from the baseline time effects of treatment III and IV, respectively, 
p-values of testing other combinations on equalities of 5i(-)i52( - )>53(') and 
<?4(-) are also reported in Table 9. 

This data set has been analyzed by Fitzmaurice, Laird and Ware (2004) 
using a random effects model that applied the Restricted Maximum Likeli- 
hood (REML) method. They conducted a two sample comparison test via 
parameters in the model for the difference between the dual therapy (treat- 
ment I— III) versus triple therapy (treatment VI) without considering the 
missing values. More specifically, they denoted Group = 1 if subject in the 
triple therapy treatment and Group = if subject in the dual therapy treat- 
ment, and the linear mixed effect was 

E(Y\b) = p x + fat + /3 3 (* - 16)+ + /3 4 Group x t 

+ /? 5 Group X (t - 16)+ + h + b 2 t + 63 (t - 16)+, 

where b = (pi, 62, 63) are random effects. They tested Hq : 04 = = 0. This 
is equivalent to test the null hypothesis of no treatment group difference 
in the changes in log CD4 counts between therapy and dual treatments. 
Both Wald test and likelihood ratio test rejected the null hypothesis, indi- 
cating the difference between dual and triple therapy in the change of log 
CD4 counts. Their results are consistent with the result we illustrated in 
Table 9. 
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APPENDIX: TECHNICAL ASSUMPTIONS 

We provides the conditions used for Theorems 1-5 and some remark in 
this section. The proofs for Theorems 1, 2, 3 and 5 are contained in the 
supplement article [Chen and Zhong (2010)]. The proof for Theorem 4 is 
largely similar to that of Theorem 1 and is omitted. 

The following assumptions are made in the paper: 

Al. Let S(9j) be the score function of the partial likelihood CBj{0j) for a 
q-dimensional parameter 6j defined in (2.5), and 6jo is in the interior 
of compact Qj. We assume E{S(9j)} / if 0j / 9jo, V&r(S(6jo)) is 

finite and positive definite, and E( 9 ^g j °' ) exists and is invertible. The 
missing propensity ~Kji m {8jo) > bo > for all j,i,m. 
A2. (i) The kernel function K is a symmetric probability density which 
is differentiable of Lipschitz order 1 on its support [—1,1]. The 

bandwidths satisfy rijh 2 - / 'log 2 rij — > oo, n'J^tij — > and hj — > as 
n j — > oo . 

(ii) For each treatment j (j = l,...,fc), the design points {tji m } are 
thought to be independent and identically distributed from a super- 
population with density fj(t). There exist constants 6/ and b u such 
that < bi < sup tgS fj(t) <b u < oo. 
(hi) For each hj and Tj, j = 1, . . . , k, there exist finite positive constants 
otj, bj and T such that otjTj = T and bjhj = h for some h as h — > 0. 
Let n = Si=i n j> n j/ n — ^ Pj for some nonzero pj as n — > oo such 
that Ei=iPj = L 

A3. The residuals {£ji} and {uj{\ are independent of each other and each 

of fe, } ) and are mutually independent among different i or i, 

resnectivelv max ll/ll - a S r,^ '^)) ( ' ^ 

respectively, maxi<j< nj ||Ujj m || — o p \n- [iogrij) j, 

maxi<j< nj E\£ji m \ i+r < oo, for some r > 0; and assume that 
lim (njTj)- 1 Y, E E{X Jim X T jim } = E x > 0, 

rij— >oo * — rf ^ — * J 

i=l m=l 

where Xjj m = Xjj m — E(Kji m \tji m ). 
A4. The functions gjo(t) and /i,* (t) are, respectively, one-dimensional and 
p-dimensional smooth functions with continuously second derivatives 
on S= [0,1]. 

Remark. Condition Al are the regular conditions for the consistency 
of the binary MLE for the parameters in the missing propensity. Condi- 
tion A2(i) are the usual conditions for the kernel and bandwidths in non- 
parametric curve estimation. Note that the optimal rate for the bandwidth 

— 1/5 

hj = 0{n- ) satisfies A2(i). The requirement of design points {tjim} in 
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A2(ii) is a common assumption similar to the ones in Miiller (1987). Con- 
dition A2(iii) is a mild assumption on the relationship between bandwidths 
and sample sizes among different samples. In A3, we do not require the 
residuals {£ji} and {uji} being, respectively, identically distributed for each 
fixed j. This allows extra heterogeneity among individuals for a treatment. 
The positive definite of T, x in condition A3 is used to identify the "param- 
eters" (Pjo,'JjO:9jo) uniquely, which is a generalization of the identification 
condition used in Hardle, Liang and Gao (2000) to longitudinal data. This 
condition can be checked empirically by constructing consistent estimate 
of Ti x . 

Acknowledgments. The authors thank the referees, Associate Editors 
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sentation of the paper. 

SUPPLEMENTARY MATERIAL 

Supplement to "ANOVA for Longitudinal Data with Missing Values" 

(DOI: 10.1214/10- AOS824SUPP; .pdf). This supplement material provides 
technical proofs to the asymptotic distributions of the empirical likelihood 
ANOVA test statistics for comparing the treatment effects with respect to 
covariates given in Theorems 1 and 2, the asymptotic normality of the em- 
pirical likelihood ratio based ANOVA test statistic for comparing the non- 
parametric time effect functions given in Theorem 3 and justifies the usage 
of the proposed bootstrap procedure. 
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